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Motivated by the fundamental question of the fate of interacting bosons in flat bands, we consider 
a two-dimensional Bose gas at zero temperature with an underlying quartic single-particle dispersion 
in one spatial direction. This type of band structure can be realized using the NIST scheme of spin- 
orbit coupling [Y.-J. Lin, K. Jimenez-Garcia, and 1. B. Spielman, Nature 471 , 83 (2011)], in the 
regime where the lower band dispersion has the form Eh ~ fci^/4-|-fcy-|-..or using the shaken lattice 
scheme of Parker et al. [C. V. Parker, L.-C. Ha and C. Chin, Nature Physics 9, 769 (2013)]. We 
numerically compare the ground state energies of the mean-held Bose-Einstein condensate (BEC) 
and various trial wave-functions, where bosons avoid each other at short distances. We discover 
that, at low densities, several types of strongly correlated states have an energy per particle (e), 
which scales with density (n) as e ~ in contrast to e ~ n for the weakly interacting Bose gas. 

These competing states include a Wigner crystal, quasi-condensates described in terms of properly 
symmetrized fermionic states, and variational wave-functions of Jastrow type. We hnd that one 
of the latter has the lowest energy among the states we consider. This Jastrow-type state has 
a strongly reduced, but hnite condensate fraction, and true off-diagonal long range order, which 
suggests that the ground state of interacting bosons with quartic dispersion is a strongly-correlated 
condensate reminiscent of superhuid Helium-4. Our results show that even for weakly-interacting 
bosons in higher dimensions, one can explore the crossover from a weakly-coupled BEC to a strongly- 
correlated condensate by simply tuning the single particle dispersion or density. 


I. INTRODUCTION 

One of the most remarkable advances in ultra-cold 
atomic gases in the recent years has been the ability to 
engineer at will, dispersions with single-particle degen¬ 
eracies or almost completely flat bands. For example, 
optical superlattices have been used to generate honey¬ 
comb, Kagome and Lieb lattice geometries mm- Lat¬ 
tice shaking mm and Raman assisted tunneling in real 
and spin space has been used to realize spin-orbit cou¬ 
pling (SOC) [316], synthetic vector potentials, and sub¬ 
sequently topological bands m [HI- Attention has now 
turned to studying the interplay between these non¬ 
trivial single-particle band structures, spin and interac¬ 
tions, which paves the way to accessing a rich variety 
of phases such as skyrmion lattices |3, integer and frac¬ 
tional Chern insulators [Mill], Wigner crystals [H] , and 
other exotic states m- Here we present a variational 
study of a low-density, 2D Bose gas at zero tempera¬ 
ture in a dispersion which is quartic in one direction, 
and which can be realized experimentally 1113 . 

An interesting example of non-trivial interplay be¬ 
tween single-particle degeneracies and interactions is a 
2D Rashba SOC gas. Here the low-energy dispersion 
has an infinite ring degeneracy in momentum space, and 
the the density of states has the form dn/dE ~ , 

typical for ID systems. At low densities, atoms sam¬ 
ple the ring degeneracy and interesting physics emerges. 
The consequences of this were first explored by Berg, 
Rudner and Kivelson m in the context of Fermi gases. 
They observed that while the kinetic energy delocalizes 
the particles over the Rashba ring, atoms can minimize 
the short range interaction energy by localizing in mo¬ 


mentum space. This competition produces a plethora of 
possible symmetry broken ground-state phases ranging 
from Wigner crystals to ferromagnetic nematic states. 

Even more interesting and perhaps less understood is 
the fate of bosons in single-particle degeneracies. On the 
one hand, by developing fermionic correlations, bosons 
can completely avoid (spinless case) or suppress (spin¬ 
ful case) short-range repulsive interactions, but such a 
state is spread out in momentum space. The kinetic en¬ 
ergy cost associated with this spreading is parametrically 
lower in flat bands, and one can expect a regime of densi¬ 
ties where fermionized wave funetions have lower energy 
than a mean-field condensate. The key theoretical chal¬ 
lenge in addressing this question is that single particle 
degeneracies enhance fluctuation effects, rendering mean- 
field theory invalid, and Quantum Monte Carlo usually 
suffers from a sign problem, and can only study small 
system sizes. Progress has to be made either by guess¬ 
ing trial wave-functions or using field theoretical meth¬ 
ods which capture the low energy dynamies. For Rashba 
SOC [m and moat bands, Sedrakyan et al. [13 US] have 
proposed a composite-fermion description, which spon¬ 
taneously breaks time reversal and parity symmetry and 
has lower energy than the weak-coupling BEC. Spinless 
bosons in quartic bands of the form ~ fe"^, were studied 
recently using field theoretic techniques by the authors of 
Refs. [IZlIIHI: who proposed that condensation is strongly 
suppressed in favor of a liquid with algebraically decaying 
spatial correlations. 

Motivated by experiments, we address the question of 
fermionization versus Bose condensation in a 2D Bose gas 
in the NIST SOC [3 or Chicago shaken lattice scheme 
[1], where the dispersion can be tuned to take the form 
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~ k‘^/4, + ky + .... We compare the energy of the 
mean-field Bose condensate to several trial many-body 
states, summarized in Table |Tj (i) a Wigner crystal (ii) 
the absolute value and the square of the Fermi-sea wave- 
function (iii) the absolute value and the absolute value 
squared of the v = 1 Laughlin state (proposed by the 
authors of Refs. [mus]) and (iv) the Jastrow ansatz m- 
While all the wave-functions (i)-(iv) have an energy per 
particle which scales as e ~ in the low-density limit 
(to be precisely defined below), and are thus energetically 
favorable over the mean-field condensate (e ^ n), we find 
that the trial wave-function with the lowest energy is of 
Jastrow type, and has finite condensate fraction and true 
long-range order. 


II. THE MODEL 


We study a two-dimensional (pseudo)spin-l/2 Bose 
system with spin-orbit coupling, which was experimen¬ 
tally realized at NIST [^: 


^soc{^) 
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where /cl is the Raman laser wave-vector, flu is the Ra¬ 
man coupling strength, and cr^ y z are Pauli matrices. 
The spectrum of the Hamiltonian has two bands, for 
< iEfj, where Er = k?‘k\/2m is the laser recoil en¬ 
ergy, the lower band has two degenerate minima, while 
for > A:Eii it has a single minimum at fc = 0 [5]. 
While the dispersion around each minimum is parabolic, 
at 0/j = the dispersion in the a;-direction develops 
a quartic structure. In the case of a Bose gas, this gives 
rise to interesting behavior at low densities, which is the 
main topic of the paper. From now on, we will be in¬ 
terested only in the Via = ^En case. We remark that 
while the Rabi coupling term explicitly breaks physical 
time-reversal symmetry, this Hamiltonian has an addi¬ 
tional Z 2 symmetry associated with the transformation 

I j T) ^ I fcs I 4 -) ■ 

Expressing energy in units of Er and momentum 
(length) in units of hk^ (l/fc^), the dimensionless single¬ 
particle Hamiltonian reads: 


Hk — k^ 2kx<Jz + 2 {^cFx -f 1). (2) 


We choose the energy offset such that the minimum of 
the lower band is at zero energy. 

We assume the interactions between particles are de¬ 
scribed by a spin-independent contact potential (in units 
of Er and l/A:/,) 


Eint(r’i -r 2 ) = g <5(ri - r2) 1^10^2, (3) 


in the space of two spins, 

where Sj S (f, i). In reality, the interactions are typi¬ 
cally spin-dependent, however our results are insensitive 
to spin dependence. We emphasize that throughout, we 
focus on the regime of weak interactions, but nonetheless 
find interesting ground states by engineering the single¬ 
particle dispersion. 

The spectrum of Elk is e± = k"^ ± 2-\/fc|~TT -|- 2 and 
the lower-band energy can be expanded around fe = 0 as 
e-{k) = k^/A + ky + .... The lower-band eigenstates of 
Hk are 


's^{k) 

Kf, 

kx — •/! -f fc^ 


V k 
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where A4 = 1 + (fc^ — a/I + fc^)^ is the normal¬ 

ization factor. Notice that at low densities (n <C k\, 
in original units), particles occupy only the states close 
to the minimum of the band, i.e. the width of the mo¬ 
mentum distribution A/c^, —>■ 0 as n —>■ 0. In that case, 
Eq. Q reduces to [s^(fe) s^(fc)] = [—I \]/\/2, and spin 
eigenstates become (approximately) momentum indepen¬ 
dent. The gas then becomes effectively spinless, and is 
described by the Hamiltonian: 





(5) 


with Vint(»’i — ^ 2 ) = g 5{rI — r 2 ). Such a Hamiltonian 
can be directly realized using the shaking lattice scheme 
of Parker et al. [4] . 

In the first part of this paper, we focus on the physics of 
the effective Hamiltonian Eq. ([^ above, and then show 
that our conclusions remain unchanged even after the 
inclusion of spin (corresponding to the NIST scheme [5]). 


III. BOGOLIUBOV MEAN-FIELD THEORY 

We start by considering the most conventional descrip¬ 
tion of a 2D Bose gas at zero temperature, namely the 
Bogoliubov mean-field description. The main assump¬ 
tion in Bogoliubov’s approach is that the majority of 
particles are condensed in fc = 0 state, and others oc¬ 
cupy fc / 0 states in the vicinity. Repulsive interactions 
deplete the condensate [21], and at the mean-field level, 
the energy per particle is given by e = 5n/2, where n is 
the density. The condensate depletion is readily found to 
be [22]: 
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where g = 2mUo/ti'^ {Uq > 0 is the contact interaction 
strength and in a quasi-2D regime Uq = a/{mUz) 

where a is a 3D scattering length and Gz is the confine¬ 
ment length in 2; direction [IS]), is a unit operator 


where ek is a single-particle dispersion, riex is the density 
of depleted particles, uq is the condensate density, and g 
is the interaction strength. The behavior of the integral 
in ([^ is usually a good indication of the fate of a BEG: for 
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zero temperature, 2D and 3D systems with a parabolic 
dispersion, the integral is convergent, fluctuations do not 
destroy long-range order. In ID, it diverges, signaling 
the absence of true long-range order. 

In our case, Sk = k'^/'i + ky, the integral is convergent: 
Uex = 3 . 854 ( 5710 )^/'^ (dimensionless variables). However, 
the ratio of the number of excited and condensed particles 
JT-ex/jT-o ~ n-Q (in a usual 2D parabolic case Uex ~ ng) 
shows that the Bogoliubov approach breaks down at low 
densities. This suggests that in the low-density limit, the 
ground state is qualitatively different from a mean-field 
condensate. 

In the 3D case, Sk = k'^/A + ky + k"^, the ratio of the 
number of excited and condensed particles is nex/ng 
v^g^■ Therefore, in the low-density limit, the Bogoliubov 
description is valid, and we expect a mean-field BEC to 
provide a good description of the ground state. 


IV. WIGNER CRYSTAL STATE 

The first example of a strongly-correlated bosonic state 
we consider is the Wigner crystal (WC), proposed by 
Berg et al. m. for the Rashba SOC case. The state is 
constructed by dividing the volume (area) in an array of 
identical rectangular boxes of size Ly, and putting 
each particle in a different box. In contrast to a mean- 
field BEC state, the interaction energy of WC is zero, as 
particles completely avoid one another. This comes at the 
cost of higher kinetic energy, as single-particle states are 
localized in boxes, compared to a BEC, where occupied 
states extend throughout the entire volume. 

To calculate the energy per particle of the WC state, 
it suffices to solve for the ground state of a single particle 
in a box. The calculation for the case of quartic disper¬ 
sion is shown in Appendix A, and for Hamiltonian (§, 
it gives Eg{Lx, Ly) = 1.2Sb{'K /-|- {■KjLy'f' , where 
and Ly are the length and the width of the box. It is 
clear that at low densities (small fe), the kinetic energy 
is “cheaper” in the x than in the y direction. This means 
that we can lower the energy by deforming the box such 
that it is shorter in x and longer in y direction, while 
keeping the total volume of the box (V = L^Ly), and 
the density (n = 1/E) fixed. We find the ratio LyjLx 
which minimizes Eg{Lx, Ly) is Ly/L^ = 0.340 and 

the ground-state energy per particle is Eg = 43.5 n^/^. 
Indeed, for the spinless case, the WC has lower energy 
than a mean-field BEC at low densities. By numerically 
solving the corresponding spinful problem [Eq.([^], we 
have checked that the energy per particle is identical to 
the spinless case in the large L^ limit. 

The WC state obviously has lower energy than a mean- 
field BEC at low densities, however it is a crystalline 
state which breaks translational symmetry. While this 
is expected to happen in low-density systems with long- 
range interactions, contact interactions typically do not 
favor formation of a crystal [23]. Therefore we expect a 


strongly-correlated state which is translationally invari¬ 
ant to have even lower energy than WC state. 

We notice that here we considered only a particular 
type of a WC (rectangular lattice) and that different 
types of WC, e.g. triangular-lattice crystal, could have 
lower energy. Still, we later show that the Jastrow-type 
state has energy e = 6.6 which is smaller than our 
WC state by a factor of 7, and we do not believe different 
types of WC can achieve such low energies. 


V. STRONGLY-CORRELATED GAS IN THE 
LOWER BAND 

A. Non-interacting Fermi gas 

The Wigner crystal example motivates us to look for 
other strongly correlated states, constructed out of low¬ 
est band wave-functions. One natural way to build cor¬ 
relations is to write down wave-functions where bosons 
avoid one another at short distances. To see how this 
lowers the energy, consider a non-interacting Fermi gas 
in the single-particle dispersion of Eq. ([^ . The density of 
states corresponding to ^ is dn/dE = (3/2)/(27r)^/^ x 
r(5/4)/r(7/4) E-^/^, where r{x) is the gamma function. 
The energy per particle in the non-interacting Fermi gas 
is then e = 6.84 which is indeed lower than a mean- 
field BEC (e = 577 / 2 ), and the Wigner crystal at low 
densities. 

It is well known that in a low-density ID system with 
contact interactions, when the contact interactions dom¬ 
inate the kinetic energy, the Fermi gas has lower energy 
than the mean-field BEC at the same density. This leads 
to “fermionization” of bosons, and the formation of a 
Tonks-Girardeau gas [24] . 

We now compute the ground state energy of several 
appropriately symmetrized fermionic wave-functions. We 
first consider the spinless case, and then generalize our 
results to include spin. 


B. Spinless system 

1. “Fermionized” many-body states 

The ground state of a non-interacting Fermi gas in the 
Hamiltonian ([^ has the following momentum distribu¬ 
tion widths: Akx = (4Ef)^/^ 77 ^/^, Aky = Ep'^ ^ 

77 ^/^, where Ep is the Fermi energy. This means that, 
at low densities, the energy is minimized by broadening 
the distribution in the direction where kinetic energy is 
“cheap” (x direction) and squeezing it in the direction 
where energy is expensive (5 direction). The WC state 
discussed above has the same property: Akx ~ ^/Lx ^ 
and Aky ~ Ly ^ . 

To construct more general strongly-correlated bosonic 
wave-functions for the spinless gas, we take a fermionic 
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state with the property Akx Aky ~ and 

construct corresponding Bose wave-functions: for ex¬ 
ample ijjB = li’pl, V'-B = V'l or = IV'Fp- This 
way we obtain a symmetric bosonic wave-function which 
obeys Afc^; ~ Aky ~ and has kinetic energy 

E\iin/N ~ {Akx)'^/4, + {AkyY ~ while the interac¬ 
tion energy is identically zero by construction. (In gen¬ 
eral, we can consider higher powers ipB = IV'fI" for n > 2 
or tjjjP, for n > 1 , but these have higher energy, as dis¬ 
cussed below.) 

The total energy is then simply given by 

.^kin — j' 

where Sk = k^/^ + ky and Uk is the momentum distribu¬ 
tion, normalized so that J dk nk = N. We compute the 
momentum distribution, and calculate the energy using 
Monte Carlo integration (see Appendix B for details). 

But which fermionic wave-functions should we choose? 
A natural choice is the ground state of a non-interacting 
Fermi gas {tpp Q). ijjpfi is a real function of spatial coordi¬ 
nates, and we construct two bosonic trial wave-functions: 
V'B.abs = IV'F.ol and V’B.sq = (V'F.o)^- 

In the case of the wave-function V's.abs, the integral 
Q diverges. The reason is that the first derivative 
of the wave-function is not continuous at points where 
V'B.abs = 0, which leads to a ~ decay of the mo¬ 

mentum distribution for large |fc|. By contrast, '0 b, sq has 
a continuous first derivative, and its momentum distribu¬ 
tion vanishes for \kx\ > 2,kF,x-i \ky\ > 2kF^y (kp^x and kp^y 
are Fermi momenta of '0 f,o ia x and y direction). The 
corresponding energy per particle is e = 13.1 which 
is considerably lower than the WC energy. 

A more exotic choice is the composite-fermion wave- 
function considered in Ref. m in the context of a 2D 
Bose gas with Rashba SOC: 

0B,cf = A? I 2 ;, - Zj l exp(- \zj\'^/A), ( 8 ) 

i<j 3 

where Zj = Xjjax + iyjjay, Xj and yj are particle co¬ 
ordinates, ax and ay are length-scales in x and y direc¬ 
tion, and Af is the normalization factor. This state has 
been shown to be a quasi-condensate with algebraically 
decaying correlations [251 US] j but it does not break time- 
reversal symmetry. In order to make wave-function have 
Akx ~ Aky ~ the lengths have to scale with 
density as Ux ay ~ Once again, the 

first derivative of ipB.ci is not continuous at points where 
V'B.cf = 0, and this leads to ~ algebraic decay of 

the momentum distribution for large |fc| rendering the 
integral Q divergent. 

However, the square of ipB.ci wave-function: 
i’Zi = (V'B.cf)^ = Af'H 1^* - (9) 

i<j 

is free from these problems. It is analytic and has an 
exponentially decaying momentum distribution for large 


|fc|. Subsequently, the integral ([7) is convergent. We find 
that the choice of length scales which minimizes the total 
energy per particle is Gx = 0.55 ay = 0.29 

and the energy is e = 9.2 which is the lowest energy 
of all the wave-functions considered so far. This state has 
zero condensate fraction, and is therefore not a true Bose 
condensate. However it has algebraically decaying corre¬ 
lations p(r, r') = p{\r — r'|) = {ij;^r)^{r')) ~ I/|r — r'\ 
as |r — r'l — 00 , and is thus a quasi-condensate [251125] . 
We reproduced this result using our Monte Carlo ap¬ 
proach (see Fig. I). While it is certainly possible to con- 



FIG. 1: We show the corelation function p{r) for two differ¬ 
ent states: (a) composite-fermion state (red dash-dotted 
line) and (b) Jastrow wave-function 0j (black line). Here we 
set y = 0 and concentrate at the dependence on x. pAjj ct 
wave-function has algebraically decaying correlations, i.e. it 
has a quasi-long range order, while ipj has a true long range 
order (see text for details). 

sider even higher powers of these wave-functions 

have higher energy as the increasing exponents broaden 
the momentum distribution of the state. 

Note that even though we casually refer to the states 
V’B.abs, '0B,sq, '0B,cf, “fermionized,” the 

issue of fermionization is a subtle one. Strictly speak¬ 
ing, our ability to express a bosonic ground state wave- 
function in terms of properly symmetrized fermionic 
wave-functions does not necessarily imply that low- 
energy excitations of this state have fermionic statistics. 
To elucidate the nature of a bosonic state in two dimen¬ 
sions written in terms of fermionic fields (which can al¬ 
ways be done even for trivial ground states), one has to 
consider a gauge theory, e.g., either arising from a par- 
ton construction or Chern-Simons flux attachment (such 
as implemented by Sedrakyan et al. for the bosonic 
Rashba model). On the other hand, there usually exists 
no simple way to write the corresponding many-body 
wave-function, which would faithfully describe gauge 
fluctuations, and those may have important and quali¬ 
tative effects on conclusions of a naive mean-field theory. 
For example, the many-body wave-function ipBpf is a 
natural mean-field description of a “fermionized” state, 
where fermions, obtained from original bosons via Chern- 
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Simons flux attachment, form the integer ly = 1 quantum 
Hall state. As discussed above, the symmetrized bosonic 
wave-function, ipB,cf, does not have a long-range order 
and hence appears to describe a strongly-correlated liquid 
state with algebraic correlations or equivalently a quasi¬ 
condensate. However, the (more general) Chern-Simons 
gauge-theory of the “fermionized” state yields a differ¬ 
ent conclusion m- integrating out fermions produces 
another Chern-Simons term, which exactly cancels the 
term associated with statistical transmutation, and what 
is left is a gapless (Maxwell) theory. It corresponds to 
a Goldstone mode and indicates broken symmetry, or in 
other words a true condensate with long-range order. In 
fact, the state proposed by Sedrakyan et al. m belongs 
to this category and is a strongly-correlated BEC, rather 
than an exotic Bose liquid. 

All in all, the field-theoretical approach based on true 
fermionization of bosonic fields and the variational ap¬ 
proach involving “fermionized” wave-functions are not 
equivalent. The former provides more insight into the 
nature of excitations, but does not easily allow for a quan¬ 
titative analysis. On the contrary, the latter can be used 
for explicit calculations of energy and other observables, 
but it does not easily elucidate the nature of low-energy 
excitations. One strategy here is to start with the varia¬ 
tional approach and explore field-theoretical description, 
if any, of a “fermionized” mean-field state, if such indeed 
comes out as the lowest-energy trial state for a given 
Hamiltonian. This however does not seem to happen in 
our case, as discussed below. 

2. Jastrow Ansatz for a strongly-correlated BEC - the 
winner 

One advantage of using “fermionized” wave-functions 
to approximately describe a ground state of interacting 
bosons is that they immediately minimize the interaction 
energy for any contact interaction (for spinless bosons), 
by the virtue of the simple fact that two fermions can 
not occur in the same point. However, there exist in¬ 
finitely many wave-functions that accomplish the same, 
without relying on any fermionic analogy. Related con¬ 
structions have been discussed in the literature, notably 
in the context of strongly-correlated BEC in Helium-4. 
Inspired by these previous studies, we now consider a 
Jastrow ansatz [28] of the following form: 

N 

V'J =A/'"I|<^(n-rj), (10) 

i<j 

^{r) = 

where Af" is the normalization, and bx, by are parameters 
describing correlation length-scale in x and y direction. 
The density, n = N/V is another important parameter 
of wave-function (jlOj). Jastrow-type wave-functions are 


generally very good at capturing the behavior of Bose 
gases ranging from small to large scattering lengths, i.e. 
from a weakly interacting to unitary regime |19j . A key 
difference here is that while usually, the Jastrow form is 
used to capture the short distance structure of the two- 
body wave-function on length scales comparable to the 
true atomic potential, here we work in a regime where bx 
and by are on the order of the inter particle spacing, thus 
much larger than the scattering length. Our ansatz is 
therefore phenomenological in nature, and does not stem 
from a microscopic calculation of the two-body problem. 

As with the previously considered trial states, the Jas¬ 
trow wave-function has the property that its interaction 
energy is zero. By choosing bx by ~ we 

can “squeeze” the system in the x and “stretch” it in the 
y direction, so that e ^ We find the optimal param¬ 
eter values are bx = 0.66 by = 0.29 and the 

energy is e = 6.6 The Jastrow wave-function there¬ 

fore has even lower energy than the composite-fermion 
wave-function V'b cf plot the single-particle 

density matrix p{r,r') as a function of \r — r'\ corre¬ 
sponding to and the Jastrow wave-function found 

above. Indeed the Jastrow form has true long range or¬ 
der, and describes a Bose-condensate with condensate 
fraction ng/n = 0.74. 

We therefore conclude that for the spinless Hamilto¬ 
nian [Eq.([^], although bosons can lower their energy by 
developing short range correlations, the correlations are 
not strong enough to completely destroy BEC at zero 
temperature. 

The ansatz wave-functions that we considered all have 
the property that ip = Q when = Vj which means 
Aint = 0. While this should be true in the n —> 0 
limit, at finite densities we expect the interaction en¬ 
ergy not to be strictly zero. In Appendix B, we estimate 
that for small densities Eint/N ^ which means 

that Eint/Ekin Therefore, as in the Lieb-Liniger 

gas [55], at low densities ^ E]^in- 

C. Spinful system 

We now turn our attention to the spinful Hamilto¬ 
nian 0 which corresponds to the NIST SOC scheme, 
and ask whether our conclusions remain valid in this case. 

We start by writing the spinful state \ipB,s) as: 

IV'B.s) = ^ /^(fcl, feAr)|fel---fcy)s, (11) 

ki. ..kN 

where 

fB{ki,...,kN) J dri...drN il;B{ri,...,rjv) 

X g-iiki-ri + ...+kN-rN) 

is the Eourier transform of the spinless wave-functions 
ipB considered above. Here |fci...fcy)s = |fei)s C ... C 
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|fcAr)s, where |fc)s is a lower-band eigenstate of Q. We 
therefore construct a spinful state exclusively from lower- 
band eigenstates. A similar construction was applied in 
Ref. [TS]. 

At low densities, the lower-band spectrum of the spin¬ 
ful Hamiltonian is the same as the spinless dispersion ([^. 
Since, by construction, the spinful state [Eq. 0 ] has the 
same momentum distribution as the corresponding spin¬ 
less state, their kinetic energy is the same (see Appendix 
D). 

However, the more complicated question is: what is 
the interaction energy of the spinful state? Since we ex- 
plicily construct the spinful many-body state only from 
the lower-band single-particle states, it is impossible to 
satisfy = Vj,...) = 0, V(j,j), for all the differ¬ 
ent spin components ■..,tw) [(Jj G 

Therefore, unlike in the spinless wave-functions consid¬ 
ered previously, the interaction energy will be finite. Still, 
in the low-density limit, we expect the zero overlap con¬ 
dition ('0 = 0 when = Vj) to be almost satisfied, since 
the system is almost completely polarized. We thus ex¬ 
pect the spinful state to have a very low interaction en¬ 
ergy. 

The interaction Hamiltonian (|^ is diagonal in real 
space, and to calculate the interaction energy it is use¬ 
ful to find a real-space representation of |'0 b,s). Unfor¬ 
tunately, the real-space representation is quite cumber¬ 
some: there are 2^ spin components (although only A^-l-l 
of them are independent due to the symmetric nature of 
/b), and expressions are difficult to obtain: 

fci.....fciv (13) 

where Sa^{k) are given in Q. 

In Appendix E we present the method to estimate the 
interaction energy, and we show that for wave-functions 
V'B.sqj V'^cf> Jastrow wave-function the energy is 

E'mt/N ~ at low densities. Therefore, Eint/Ekin —>■ 
0 when n —)■ 0 and Etot/N = (Ekin + Eint)/N —>■ 6.6 
for the Jastrow state. 

It is important to assess the validity of constructing 
the spinful state only from lower-band eigenstates: we 
have already shown that the ground state energy cannot 
be greater than e ~ If there was a hnite fraction u 

of particles occupying the higher band as n —>■ 0 , then the 
energy would be E/N ^ uA, where A is the gap between 
the two bands. However, this clearly contradicts the fact 
that E/N < . Therefore, m —>■ 0 as n —0 and the 

n —>■ 0 ground state will only contain states from the 
lower band. 


Wave-function 

e Section 

Mean-field BEC 

gn/2 


III 

Wigner crystal 

43.5 


IV 

Absolute value of Fermi-sea w.f. 

OO 

V 

B1 

Fermi-sea w.f. squared 

13.1 nU3 

VBl 

Composite-fermion w.f. 

OO 

VBl 

Composite-fermion w.f. squared 

9.2 

VBl 

Jastrow state 

6.6 

VB2 


TABLE I: Energy per particle (e) of different states in the 
low-density limit (n and g are dimensionless density and in¬ 
teraction strength, respectively). Of all the wave-functions 
(w.f.) we consider, the Jastrow state has the lowest energy. 
Two wave-functions (absolute value of Fermi-sea w.f. and 
composite-fermion w.f.) have diverging expectation value of 
k/. (see text for details). 


VI. DISCUSSION AND EXPERIMENTAL 
RELEVANCE 

In this paper, we considered a system of interacting 
bosons with a quartic single-particle dispersion. It was 
shown that the low-density limit of the model hosts a 
strongly-correlated ground state, where the mean-field 
Bogoliubov state can be easily ruled out as being para¬ 
metrically higher in energy than the strongly-correlated 
states, where bosons develop local correlations and avoid 
each other. 

Among the many trial states we considered, a long- 
range-ordered condensate described by the Jastrow wave- 
function [Eq. 0 ] was found to have the lowest energy 
per particle of e = 6.6 (compared to e = gn/2 for 
a mean-field BEC). This is in agreement with Ref. [17] 
where it was argued that the ground state of system © 
has long-range order. The condensate fraction was found 
to be Nq/N = 0.74. i.e, it is a strongly-correlated BEC 
with significant depletion of the condensate due to the in¬ 
terplay between interactions and the unusual band struc¬ 
ture. 

Importantly, the mean-field BEC and the Jastrow BEC 
break the same symmetry, therefore, we expect that the 
system continuously evolves from a weakly to a strongly- 
correlated BEC state from high to low-densities, with¬ 
out any phase transitions in between (a similar weak- 
to-strong-coupling cross-over can be tuned by evolving 
the single-particle dispersion from the usual quadratic to 
quartic). 

Note however that in the absence of a systematic pro¬ 
cedure to explore many-body ground states of strongly- 
correlated systems, our variational-approach results are 
strongly suggestive, but not conclusive. Eventually, it is 
experiment that would fully elucidate the nature of the 
ground state, and to realize our model is at the experi¬ 
mentalists’ fingertips. 

The strongly-correlated condensate we predict can be 
detected experimentally using a number of probes. For 
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example, the suppressed condensate fraction is measur¬ 
able in time-of-flight [30]. Another signature of strong 
correlations is the ratio of interaction and kinetic energy 
which is very small at low densities. This could be ac¬ 
cessed via quantum quench experiments, i.e. the inter¬ 
action parameter g could be suddenly changed and the 
effect on the total energy of the gas could be measured. 
Strong local correlations can also be measured in situ 
by observing the anti-bunching of bosonic atoms m- 
Finally, several groups [SS] [33] have directly measured 
p{r,r'). This would give information about the conden¬ 
sate fraction, and the type of order present in the gas. 

We can estimate the density below which strong cor¬ 
relations become energetically favourable by equating 
mean-field-BEC and Jastrow-state energies (see Table [l|. 
In the case of ®^Rb with a z-direction confinement fre¬ 
quency ujz = 2tt X 4000 Hz, this gives n « 10“® « 

6 X 10^ m“^ (/cl = 27r/A, where A « 800 nm [34]), which 
is much lower than typical densities in cold-atom exper¬ 
iments studying 2 D systems {n ~ 10^^ m“^ [Si]). How¬ 
ever, using Feshbach resonances to increase 5 , it is pos¬ 
sible to make strong correlations favourable at consider¬ 
ably higher densities, up to n « 0.004 « 2 x 10^^ m“^, 

which could be achieved experimentally. At densities 
higher than this the dispersion in x-direction cannot be 
approximated by a quartic term anymore, and higher- 
order terms have to be included. 
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Appendix A: Box-potential ground state in a system 
with quartic dispersion 

Here we show how to find a spectrum of a particle in 
box potential with Hamiltonian H 4 = = d'^. While 

H 4 is similar to the usual quadratic dispersion in a sense 
that both are diagonal in momentum space, there is one 
fundamental difference: in a system with H 4 , not only 
the wave-function, but also its first derivative has to be 
continuous for the wave-function to have finite energy 
expectation value. 

If we choose a box of length L and — L/2 < x < L/2, 


then the boundary conditions are ^/>(—L/2) = ip{L/2) = 
0 and dx'>p{—L/2) = dx4’{Lt2) = 0. Solutions of the 
equation d^ip = Eip are exp(fcx), exp(—fcx), exp(zfcx) 
and exp(—z/cx), where k = Since iJ 4 is symmetric 

under inversion (x —^ —x), we expect a symmetric ground 
state: 

■00 (a:) = oi -I- 02 cos(fcx), (Al) 

where oi, 02 are coefficients that have to be deter¬ 
mined. Boundary conditions then require tan(fcL/2) = 
— tanh(fcL/ 2 ), which can be solved graphically: in the 
ground state kL = 4.730 and Eq = 5.140 (tt/L)'^. The 
ratio of coefficients is 02/01 = 15.06. 


Appendix B: Monte-Carlo calculations 

Here we describe Monte-Carlo methods we used to cal¬ 
culate the kinetic energy of various trial wave-functions. 
We were primarily interested in finding the expectation 
value of Hamiltonian (§ and for analytic wave-functions 
this can be done in two ways. As in standard Variational 
Monte Carlo techniques, the first step is to sample the 
“local energy”, Eioc = [36] : 

~ /|0|2dJ2 
/|0|2di^ ’ 

where H = ~ i dR = dridr 2 ■ ■ ■ dr^. 

We then use a Metropolis algorithm to sample the 
local energy with probability distribution P{R) = 
|0(fi)IV/|0(fi')|2dfi'. 

In the case of non-analytic wave-functions like 0_B,abs 
and tpB.ci the expectation value of dx^ cannot be calcu¬ 
lated this way because a finite energy is associated with 
points which have discontinuous derivatives of 0. The 
correct method in that case is to first calculate the mo¬ 
mentum distribution of the state, and then compute the 
expectation value oi Sk = kx/4, + ky. We calculate the 
momentum distribution n(fe) using: 

n{k) = N y'dr 2 ---drjv|/(fe,T- 2 ,...,rAr)p, (B2) 

/(fe,r 2 ,...,rjv) = ^ J (irie"*'“'”i0(ri,...,rjv), 

where we chose the following normalization: 
/|0(ii)pdi2 = I and / n(k)dk = N. This can be 
written in the form suitable for Metropolis importance 
sampling: 

n{k)=N J dR\i;{R)\^\fM{k,r 2 ,...,rN)\^ (B3) 




/Ar(fc,r 2 ,...,rAr) = ^ J drie ...,rN), 

V>(ri,...,rjv) 
dr'\tp{r',r 2 , ...,rAr )|2 

For analytic wave-functions, both methods produce the 
same result. 

The wave-functions we considered all have the property 
of being very anisotropic, i.e. they are given in terms of 
length-scales ay where Gx ^ ay. The best way to do 
calclations is then to rescale the coordinates so that in the 
new units ax ^ ay 1. For example, in the case of V'b cf 
we first calculate expectation values = (k^) and q ;2 = 
(ky) for the wave-function with ax = ay = 1 (rescaled 
wave-function). The expectation values corresponding 
to the wave-function in the original units are then simply 
04 / 0 ^ and 02 / 0 ^, respectively. In the end, we minimize 
E{ax,ay) = {kl)/d + (k^) = aA/[^a%) + a 2 /ay, while 
keeping the density n = l/(27raa;ay) constant. 

In the case of wave-functions where we can apply peri¬ 
odic boundary condition (V's.absj '0B,sq) and V'j) we did 
calculations with TV = 400 particles. However, in the case 
of composite-fermion wave-functions {ipB,ci and V'b cf) 
did calculations with N = 1600 particles. There the den¬ 
sity of a wave-function with finite number of particles has 
a form of a droplet with radius R = ^2(TV — 1) (when 
Gx = Gy = 1) and the presence of the boundary increases 
the value of finite-size correction. Larger system sizes 
were therefore necessary. 



Appendix C: Estimating the interaction energy of a 
spinless gas at small, but finite densities 

In order to estimate the interaction energy at small 
densities, we can make a simple order-of-magnitude cal¬ 
culation: we choose some coordinates A = (r 3 ,...,rjv) 
and we keep them fixed (see Appendix E for more 
details). Now we can write the wave-function as 
'!/i(ri, r 2 ; A), and we can define parameter C as a mea¬ 
sure of a ri = r 2 wave-function amplitude: 


Appendix D: Kinetic energy of a many-body 
wave-function 


Here we show that spinful wave-function constructed in 
Eq.@ has the same kinetic energy as the corresponding 
spinless wave-function. 

We consider the following state: 


= ^ /^(fci, fcAr)|fci---feAf), 


(Dl) 


where /b is normalized: J2ki,...,kN 
\ki...kiy) = |fci) (g) ... (g) |feAr) is an orthonormal 
momentum-eigenstate basis. Here state |fc) can describe 
either a spinless or spinful [Eq.Q] single-particle mo¬ 
mentum eigenstate. The state is an eigen¬ 

state of kinetic energy operator with energy eki...kN = 
Eki + ■■■ + £kN- Therefore 

Akin ~ ^ ^ ...fciv I/b (^ 1 ; ■■■; ^v) I 

ki...kN 

= (£fci + •■• + Efew) l/B(fel) feAf)P 

= X! l/B(fci,---,fcAr)P H- 

fei k2...kN 

— ^ ^ ^k^kf 
k 

where we assumed /b is symmetric with respect to par¬ 
ticle exchange and n/. is the single-particle momentum 
distribution: 


nk = N y] |/B(fe,fc2,--,fcAf) 

/C2--feiV 


(D3) 


It is clear that kinetic energy does not depend on whether 
\'ip) [Eq.(Dl 1 ] describes a spinless or spinful state, as long 
their momentum representation fs and dispersion £k are 
the same. 


Appendix E: Estimating spinful state interaction 
energy 


c = v 


J dr\'ijj{r,r)\ 
f dridr2ltp(ri,r2)l’ 


(Cl) 


where V is the volume. We first estimate the kinetic 
energy: Ekin(C') ~ — C). The reasoning is that 

for (7 = 0, Akin When (7 = 1, the gas is not 

correlated and Akin ~ 0. Moreover, Akin should not have 
an extremum around (7 = 0, and therefore should be 
linear in (7 in that region. 

The interaction energy is Aint ~ N'^g f c?r|'0(r,r)p ~ 
gnC^. We minimize Akin + Ajnt with respect to (7 and the 
optimal (7 is (7 ^ jg^ and Aint jg. This means 

Aint/Akin that is the kinetic energy is a dominant 

part at low denities. The same reasoning gives the correct 
density scaling of Aint in the low-density regime of a ID 
Lieb-Liniger gas. 


The spinful state is defined as 

IV'B.s) = y] fB(klT..,kN)\ki...kM)s, (El) 

fcl.. .fciv 


where \ki...kM)s = |fei)s < 8 ) < 8 ) |feAr)s and \k)s is a lower- 

band single-particle state. 

The real-space representation of \ki...k]\[)s is 


(n- 

• rAr;tt •• 

t fcl...fcAr)s 


S^(fcl)s^(fe2)...S^(/Cjv) 

(n- 

•rAr;tt •• 

i\ki...kN)s 


S^(fcl)s^(fe2)---S;(/Cjv) 

(n- 

.rjv; 44- 

t fci...fcAr)s 


Si{ki)si{k2)...s^(kN) 

_{ri. 

4-4-.. 

1 

-a 

— i- 


_si{ki)si{k2)...Si(kN)_ 
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^ gi(fci-ri + ...+fejv'r’]v) (E2) 

where s^{k), s^(fc) are given in eq.Q. The real-space 
representation of the spinful wave-function is then 


,rN) 

,rN) 


= ^ fsiki, ...jkisr) 

V'U-T(^ir-- ,rN) 

- ,rN)_ 


st(fei)st(^ 2 )---St(^Af) 

s^{ki)s^{k 2 )...si{kN) 


X 


^z(fci-ri + ...+fciv-riv) 


(E3) 


si{ki)si{k 2 )...s^{kN) 

si{ki)si{k 2 )...si{kN) 


We are interested in the low-density regime and there 
|fc| ^ 1. We can then expand spin coefficients as: 

Also, we can replace with —idx- For example, the 
component is then: 



X • •• X ^1 - -f ^ -h ■ ,rN), 

(E5) 

since by definition of Jb (Eq. 0 ): 

tpBiri,-■ ■ ,rN) = fsiki,--- ,kN) 

fci---/civ (E6) 

X gd'=l'’’l-l- y-kN-rn) 


get a two-body wave-function from which it is easy to 
calculate {Vi 2 ). Later we show that almost any choice of 
A gives the same value of (^ 12 )- 
We start by defining 


idx 


92 


$(ri, r 2 ; A) = N{\) 1 - ^ + ^ + ... 


1 - 


id. 


X 2 

2 


d^ 

^X 2 


X(ri,r 2 ; A), 


(E7) 


where 


X(ri,r 2 ; A) = ( 1 - 


idx 




-b ■ 


X ••• X 1 - 


idx 


'-'xn 


+' 


^B, 

(E 8 ) 

(E9) 


and A/’(A) is such that 

j dridr2 |d>(ri, r2; A)p = 1. 

We notice that $(ri,r 2 ; A) is simply with a differ¬ 

ent normalization [see Eq. ®]. 

If is analytical, we can Taylor-expand around ri — 
r2 = 0 : 

V'B =^axx {xi - X2f + ^ayy (j/i - y2f 
+ axy (xi - X 2 ) (yi - 2 / 2 ) H-, 

where we assume ipB = 0 when ri — r 2 =0 and a^- = 
aij(fcm,A), where rcm = Vi + r 2 . x retains the same 
structure, but with different coefficients: 


X =^bxx (xi - X2f + ]^hyy (yi - y2f 

+ bxy (xi - X2) (yi - 2/2) H-, 


(Ell) 


However, once we act on y with si{—idxi)si{—idx 2 ), 
function $ will have non-zero value for ri = r 2 which 
will give rise to finite interaction energy. 

Let us make the change of variables: Xr = xi — X 2 , 
Xcm = xi + X 2 - Then: 


‘l*|ri=r 2 


[. idx, dl^ 

^ 2 8 

^ (l-^-b^) x(ri,r 2 ;A)|^i =^2 


The strategy for calculating the interacting energy 
is: (a) we concentrate on expectation value of V 12 = 
g5{ri — r 2 ) (the total interaction energy will then be 
N{N — l)/2 times that value), (b) We first calculate 
the contribution to interaction energy coming from V't'-q 
component, (c) We show that all other spin components 
give approximately the same contribution. 

Estimating contribution .— The idea is to choose 
some random values for coordinates 1 ^ 3 , • • • , r jv and keep 
them fixed [we define A = (ra,- • • ,rAr)]. This way we 


1 ~ '^^cm) + gidxr + dxcntY 


1 ~ 2 dxc,n) + gi~dxr. + dx^^y 


X x(n,'r 2 ; A)| T*i—r 


= {^+\dl^ + ---^x{ri,r 2 \>y\r^=r 2 

_ ^XX 
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We can estimate b^x ~ where Akx is the mo¬ 

mentum width in x direction and |<h| is the average mag¬ 
nitude of $: 

= ^ / (E12) 

where V is the volume. The interaction energy corre¬ 
sponding to $ is then 

Ei2($) = / dridr2 g6{ri - r2)\^{ri,r2)\‘^ 

■> _ (E13) 

^Vg{Akxf\^^ ^ ^{Akxf 

It is clear that a different choice of A would give the same 
estimate. The same is true for different spin components 
V'o'i,...,o-jv Therefore, when we average over all A and 
{cti, ..., (Ttv}: 

Pai,...,aN (A)Ei2[$ 

cri,...CTjv (A)], 

<yi,..-,<7N 

P(7i,...,crjv (A) = j dridr2\'4) (ri,r2, A)|^ 

-fei4> 
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